Best linear forecast of volatility in financial time series 
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£h ' The autocorrelation function of volatility in financial time series is fitted well by a superposition of 

several exponents. Such a case admits an explicit analytical solution of the problem of constructing 
the best linear forecast of a stationary stochastic process. We describe and apply the proposed 
analytical method for forecasting volatility. The leverage effect and volatility clustering are taken 
into account. Parameters of the predictor function are determined numerically for the Dow Jones 30 
Industrial Average. Connection of the proposed method to the popular ARCH models is discussed. 

^ ' I'ACS numixTs: N<ur>.( ;h. *•!. ".].>.>. 02 .-">(). N - 02. .">(). Sk 



I. INTRODUCTION 



The methods developed in studying complex physical systems have been successfully applied throughout decades 
to analyze financial data Q,Q,IS- The quantitative study of financial data continue to attract the growing interest 
motivated by the existence of universal features in the dynamics of different markets, such as power-law tails of 
the return distributions 0, 03, IE Is l W i [lij . [TH . scaling as a first approximation Q and deviations from scaling 
of the empirical return distributions [a. \8L Jl2 | . volatility clustering J l2. Il3j . and leverage effect [Til Il5l Hi} . The 
phenomenological and microscopic models llll 13 HE EH E3 have been proposed to explain the established 

stylized facts. The field of research connected to modeling financial markets has been named Econophysics. 

A stock's volatility represents the simplest measure of its riskiness or uncertainty. Formally, the volatility is the 
annualized standard deviation of the stock's returns during the period of interest. The random walk model proposed 
by Bachelier in 1900 year p] presupposes a constant volatility. There is an ample empirical evidence, however, that 
the volatility is not a constant, but represents a random variable. Two well established stylized facts concerning 
the volatility are long ranged volatility-volatility correlations that are also known as volatility clustering 0] an d 
return- volatility correlations that are also known as leverage effect ES • 

The volatility is a key variable to control risk measures associated with the dynamics of prices of financial assets. 
The implied volatility extracted from options prices represents a market estimate of future volatility. A pure exposure 
to future volatility is provided by the volatility swaps The volatility enters all options pricing models, so its 

knowledge has a great value for estimate of the equilibrium options state-price distributions. 

The volatility clustering manifests itself in the occurrence of large changes of the index at neighboring times 
(observed localized outbursts). The leverage effect has its origin in the observed negative correlation between the 
past returns and future volatility. The possible explanation to this effect 0, 0, is due to the fact that negative 
returns increase financial leverage and extend the risk for investors and thereby a stock's volatility. A statistical study 
|23j demonstrates clearly that the leverage effect is one-directional: past returns correlate with future volatility only. 

In this paper, we propose an analytical method to evaluate future volatility as a linear function of the lagged 
volatility and lagged returns. The method takes the volatility clustering and leverage effect into account and provides 
O ■ for stationary stochastic processes the smallest forecasting error in the class of all linear functions. In this precise 
\ sense, we talk on the best linear forecast (BLF) of the volatility, 
r ■ The BLF problem for a stationary stochastic process was formulated and solved by Kolmogorov [24| in 1941 year 
and Wiener j^j in 1949 year. A modern review of the BLF methods can be found in Ref. j2(|. We apply these 
methods to construct the BLF volatility function for the Dow Jones 30 Industrial Average (DJIA). 

The outline of the paper is as follows: In the next Sect., we remove the leverage effect from the original time series 
to work with a reduced volatility x(t) that has by definition a vanishing covariance with the past returns. The spectral 
density of a stochastic process can be factorized, f(u>) = \(p(uj)\ 2 , if its correlation function represents a superposition 
of the exponential functions. An explicit expression is derived for the amplitude (p(w). The analytical properties of the 
amplitude tp(u>) in the complex w-plane are important to provide an explicit representation of the predictor function. 
In Sect. 3, the BLF problem is analyzed further to account for the reduced volatility clustering and to construct the 
BLF function. In Sect. 4, we fit 100+ years of data of the daily historical volatility of the DJIA in order to determine 
parameters of the BLF function. Numerical estimates are given to illustrate the developed method. The minimization 
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of the forecasting error for the reduced volatility predictor function is shown to be equivalent to the minimization of 
the forecasting error of the original volatility time series. An explicit expression for the forecasting error is given. In 
Conclusion, a connection of the BLF method with the ARCH models [13, E3, US H3 > m which future variance is also 
represented as a linear combination of the past observables, is discussed. 



II. FACTORIZATION OF SPECTRAL DENSITY 



The evolution of a market index value or a stock price S(t) is described by equation (see e.g. |3J): 

dS(t) 



S(t) 



(j,dt + dip{t). (1) 



The value dip(t) is a noise added to the path followed by S(t), with the expectation value E[dV(0] — an d the variance 
of V&r[dip(t)] — a(t) 2 dt. The volatility a(t) represents a generic measure of the magnitude of market fluctuations. We 
consider a discrete version of the random walk problem by setting dt — 1, dip(t) = £(i), and dS(t) = S(t) — S(t — 1). 
The sampling intervals are enumerated by integer time parameter t. 

The volatility a(t) is a hidden variable and its extraction form the market observables is a separate difficult task. 
The possible estimator rj(t) = of the volatility is defined in terms of returns 

m = (S(t)-S(t-l))S(t)- 1 -^. (2) 

In what follows, the term "volatility" refers to the estimator r)(t) = |£(t)|, the annualizing factor will not apply. A use 
of the variance estimator v(t) = |£(i)| 2 would complexify the problem due to divergences connected to the existence 
of power-law tails ("variance of variance" is infinite, Var[£ 2 ] ~ oo, since dF(£) ~ d£/£ 4 at £ ^> 1, see e.g. 0). At 
large time scales, different estimators arc expected to be close to volatility a{t) and to each other. The problem of 
efficiency of volatility estimators is postponed for other studies. 

It is usually assumed that financial time series constitute stationary stochastic processes the autocorrelation func- 
tions of which depend on the relative time only. The stock evolution problem is assumed therefore to be invariant 
with respect to time translations. 

First, we remove from the time series rj(t) the leverage effect using the variable x(t) : 

X(t) =r ? (t)-^Cov[^(0),^(- S )]Var- 1 K]^(t- S ). (3) 

s 

The decomposition has a predictive power, since Cov[rj(t), £(t — s)] ~ 0(s), so x(t) depends on the lagged price 
increments only. Note that E[?y] = E[x], since E[£] = 0. Due to the definition and in virtue of equation 

Cov[C(*),e(s)]=^Var[e], (4) 

that holds true for sampling intervals greater than 20 Min. 0, El i we have 

Cov[ X (t),Z(s)]=0. (5) 

The reduced volatility x(i) does not experience the leverage effect. So, its predictor depends on the past only. It 
is possible therefore to focus on the volatility clustering only, while the leverage effect is taken into account explicitly 
through Eq.@. 

The autocorrelation function B(t) = Cov[x(i), %(0)]/E[x 2 ] of a stationary stochastic process x(i) can be fitted in 
many cases by a superposition of exponents 

n 

B(t) = Y / d i e- a ^. (6) 

i=l 

The best linear forecast of the observable x(t) in such a case simplifies substantially. The case n = 1 is discussed in 
Ref. |26(. We provide a solution of the BLF problem for arbitrary values of n. 

The spectral density f(w) of the stochastic process x(^) is given by the Fourier transform of the correlations 
coefficient (©: 



f(u)= V e-^ t B(t)= s Td l - £ ^ (7) 

t= — OO 1=1 v /v u> 
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where u = exp(— iui). The function f(tv) can be represented in the form 

/H = P n _^) (f[(l - e -<*«)(l - e-«4)) • (8) 



\i=l 



where (f> = h(u + —) and 



P n _i (0) = 2" exp(- y^Qjj)y^ sinh(ai) J^[(cosh(a l ) - 4>) 

7=1 7=1 k=jti 

77—1 71—1 

= A l 2"- 1 exp(-^^)n( cosh (^)-^)- (9) 

7=1 7 = 1 

The amplitude <p(u) such that f(u>) = ip(u)ip(u)* can be chosen to be analytical, rational and regular at \u\ < 1: 

71-1 / 71 \ _1 

<pm = di/ 2 n a - - e ^ u ) ■ ( io ) 

7=1 \i=l / 

The additive representation 



7=1 

is completely equivalent to the multiplicative representation (|10|l . Here, = e~ Qi and 

77 77 — 1 71 



D n = 2 exp(- «j + Vi) rfj sinh(ai), 



7=1 7=1 7=1 

77-1 




k=l 

cosh(^) = 4 (12) 

where (/)q are n — 1 roots of equation P„_i(</>q) = 0. For n = 1, Ci = c^ 2 -\/l — The analytical solutions for exist 
up to n = 5. 

The knowledge of the Fourier transform of the autocorrelation function is not sufficient for a complete reconstruction 
of the Fourier transform of the stochastic process x(t). If the autocorrelation function represents a superposition of 
the exponents I©, the problem admits a solution f{u), such that f(u) = |<£>(it)| 2 , in the class of rational functions. If 
we require further that the function ip(u) be regular at |u| < 1, an unambiguous solution ip(u) can be provided. This 
solution coincides with the Fourier transform of the time series x(t) up to a phase factor. It is remarkable, that we 
need not to know the phase, since all the relevant information is contained in the spectral density f{oj). The BLF 
problem simplifies then considerably due to the special analytical properties of the function tp(u). 



III. BLF FUNCTION 



The correlation function corresponding to the spectral density JJJ can be found from the inverse Fourier transform. 
In case of (|10pi . we consider first t > 0: 

B(t) = ^e-^)- = ^iuM-)- (13) 

where C r = {e~ zu ' ,u = — n...ir}. The poles of <p{u) are located at \u\ > R = min{e Qi }, while the poles of <p(— ) are 
located at \u\ < 1/R. We move the contour C r to infinity and get 

77 77 
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Comparison with Eq.JBJ gives 

n 

■Pip, 



k— 1 



The same result (JUJl comes out at i < 0. For n = l, B(t) = P^cf/(1 - /3f). 

A stochastic process x(^) can be represented as a linear combination of a normally distributed uncorrelatcd sequence 

T X ) With (j2 



C(t) ~N(0,a x ) with 4 = E[ X 2 ], 



+ OC 



x(t) = E[ X ] + ^C( S )C(t- S ), (16) 

s=0 

provided that the spectral function admits the factorization and the amplitude ip(u) is regular at \u\ = 1 < R (see 
e.g. 26]). The expansion coefficients equal 

n 

C{t) = Y,*Pl (17) 

i=l 

It is remarkable that only retarded £(s) enter the summation in Ea. (|16|l . This is a consequence of the convergence of 
the Taylor expansion of the amplitude f(u) at |tt| = 1, which is in turn a consequence of the analyticity at \u\ < R: 
The convergence radius of the expansion is associated with the first pole at \u%\ = R. The stationary stochastic 
process x(t) can be interpreted as a result of filtering the normal sequence ((t). 
The BLF function for the time horizon t has the form |2(j 

X T (t) = E\x] + J2^—(x(t-s)-E[ X }). (18) 

S — T 

The weight coefficients S r (0)^ are derivatives of the function S T (u) = (p T (u)/ip(u) at u = 0. Here, 

+00 



^( 11 )=^c(tK = ^ Cl (^r T i 



PiU 

For constructing a linear prognosis function, the overall normalization factor in B(t) is not important, since it drops 
out from the ratio (p T (u)/(p(u). In virtue of Eq.l|H|l. 

Cov[xr(*U(s)] = 0. (19) 

In terms of the stochastic process ("(£), the BLF function looks like 

+00 

*r(«)=E[x]+£c(*X(t-*). (20) 



At s = r, we obtain 

and at I = s — r > 

H r (0)(* 



5r(0)^ _ S?=iCie-^ r . n 



= £>- 1 /^ e - 3 a-D TJ( c -"i - e -"«) £ Ci e-^n(e-^ - e^). (22) 



S! 

j — 1 \ / i— 1 fc^i 



The last two equations complete solution of the BLF problem for the case when the correlation function is a 
superposition of the exponent functions J^J. The n(n — 1) terms in the right side of Ea. (|22|l are not all positive 
definite. 
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FIG. 1: Empirical correlation coefficients Cair[r](t), 17(0)], Corr[?7(£), £(0)] and Corr[x(t), x(0)L Corr[x(£), £(0)] versus the num- 
ber of trading days t of the Dow Jones 30 Industrial Average. The reduced volatility xiP) is defined through Eq.J3J. The 
correlation coefficients are calculated using 100+ years of the daily quotes, starting on the May 26, 1896 and ending on the 
December 31, 1999 (i.e. a total of 28507 trading days). The values of Corr[?7(i), 7/(0)] and Corr[?/(f), £(0)] are denoted, respec- 
tively, by triangles and diamonds. The values of Corr[x(i), x(0)] an d Corr[x(i), £(0)] are denoted by circles and boxes. The 
solid curve is the exponential fit © with parameters given in Table 1. It is seen that Corr[x(£), £(0)] ~ 0. The leverage effect 
is thus removed from x CI- 



TABLE I: Parameters di and a.i entering the fit of the autocorrelation function (|SJl of the reduced volatility, parameters l>i 
which determine the roots of equation P n _i(0o) = 0, and parameters d which determine the additive representation 11 U of 
the function <p(u). The value of D4 represents the normalization constant according to Eg. 1121 . 



di 0.40 


ai 


+00 v x .002257 


Cl 


.606241 


d 2 0.05 


a 2 


1/20 u 2 .012107 


C2 


.038233 


d z 0.03 


«3 


1/250 v z .125302 


C3 


.007857 


d 4 0.05 


Q 4 


1/1000 D 4 .435341 


C4 


.007473 



IV. PARAMETERS OF BLF VOLATILITY FUNCTION FOR DOW JONES 30 INDUSTRIAL AVERAGE 

Let us apply the BLF method to forecasting the volatility for the DJIA. The daily returns are defined by Eq.J2J) 
where Sit) are the DJIA index close values, the volatility equals r)(t) = |£(f)|, and x(t) 18 defined by Eq.0. In 
Fig.l, we show the empirical values of the correlation coefficients Corr[?j(i), 7/(0)], Corr[rj(t), ^(0)] and Corr{x{t), x(0)], 
Corr[x(t), £(0)] and the exponential fit of the correlation coefficient Corr[x(i), x(0)] versus the time lag t. Let us remind 
that Gorr[A,B] =VarL4, B]/^/E[A 2 ]E[S 2 ] and -1 < Cott[A,B] < 1. To calculate x{t), we run the summation over s 
in Eq.© from to 250 and use the empirical correlation function of r](t) and £(s) without additional smearing. Up to 
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TABLE II: The weight coefficients 12 U and 12211 of the BLF function for some values of the parameters / and r. Here, r is the 
forecast horizon, I = s — r is the number of trading day entering the predictor function starting from the most recent day. 



I T = l t = 2 t = 10t = 100 

.07830 .07555 .05780 .01862 

1 .06942 .06702 .05150 .01710 

2 .06158 .05949 .04594 .01575 
10 .02432 .02366 .01941 .00916 



t = 250, the correlation coefficient Corr[x(t), £(0)] is less noisy as compared with other correlators. The parameters 
di and on are listed in Table 1. The equation Pz{4>) = determines the parameters z/; (i — 1, 2, 3) according to Eq.©. 
Using Eas. l|12|l . we find the values of D4 and Cj, which we also show in Table 1. One can check that Eq. l|15fl is satisfied. 
The correlation coefficient Corr x(0)] drops from 0.53 (= J2t=i ^») to 0-13 ( = Y^t=2 ^») wnen t changes by one 
unit from t = to t = 1 . The value of a% is therefore large and can be fixed by considering high-frequency data only. 
The results shown on Fig. 1 and in Table 1 and Table 2 are obtained for ot\ = +00. The values of 1/ai for i — 2, 3, 4 
equal to about one month, one year, and four calendar years, respectively. 

The weight coefficients S r (0)( s ) can be found with the use of Eq. l|21|l at s = r and Eq.J22J) at s > r. We show 
values of the weight coefficients divided by s! in Table 2 for I = s — r = 0, 1, 2, and 10 and r — 1,2, 10, and 100. 

The BLF volatility function looks like 

+°° n (nlO) + °° 
f) T (t) = E[ X ] + ^^—^ - s ) - E W) + E Cov ^(°)' Var- x K]£(i - «) (23) 

s— r s—r 

where the unknown future returns set qual to zero: — s) — » E[£(t — s)] = for < s < t. 
Using Eqs.Q, ©, and lfTS)l. one gets 

r-l 

E[«r(t) -^)) 2 ] =E[(xr(*) -XW) 2 ] +Y, Co ^(0)^(s)?y^- 1 [^ (24) 

so the minimization of the Xr(t) error according to Ea. (|18|) is equivalent to the minimization of the fj T (t) error. Using 
decompositions (|16|l and (|20l) . the x T (t) error can be evaluated as 

r -! " n 1 _ -(a ! +Q fc )r 

E[(XrW - XW) 2 ] = E[ X 2 ] 2 C7 2 (.s) = E[ X 2 ] J2 E ^ i_ e -(a i+ a fc ) ■ ( 25 ) 

s=0 i=l k=l 

(There is a misprint in Eq.(10.2) of Ref. |2^) At r -> +00, E[(x r (*) - xW) 2 ] ^ E[ X 2 ] ^=1 rf i = Var[x], in agreement 
with the fact that X r(i) — > E[x]- The arguments of such a kind do not apply to the variance estimator v(t) = |£(i)| 2 , 
since Var[v] = 00 due to the power-law tails of the return distributions (jj. 

Nonlinear models for volatility forecasting 0, [2(| , which take into account besides the volatility clustering and 
leverage effect also heavy tails of the returns distributions and the approximate scaling, represent an alternative 
class of the stochastic volatility models. The efficiency of such models can be tested in general using Monte Carlo 
simulations and/or backtests over historical data. The approach of Refs. mil is more general, since it allows a 
calculation of the probability density function of the volatility. The BLF method predicts the average volatility only. 
It can, however, be extended to forecasting |£| a for arbitrary < a such that E[|£| 2a ] < 00. If all moments E[|£| a ] 
of the future distribution are known, the reconstruction of the probability density function of the volatility must be 
possible within the BLF method also. 



V. CONCLUSION 

The BLF problem for a stationary stochastic process was formulated in 1941 year by Kolmogorov [24[ and later by 
Wiener |25| . A modern review of the BLF methods can be found in Ref. |26| . In this paper, we reported an explicit 
analytical solution of the BLF problem for practically important case when the autocorrelation function represents 
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a superposition of exponential functions. The autocorrelation function of the volatility in a financial time series is 
known to be fitted well by such a superposition. We applied the obtained results to construct the BLF volatility 
function for the DJIA. 

The po pul ar autoregressive conditional heteroskedasticity (ARCH) models of time dependent volatility proposed 
by Engle |27j (for a review see |28l l29l l3Cjp. describe the variance a(t) 2 as a linear function of the past observables. 
The ARCH models are conceptually very close to the BLF approach. Eq. (12:>ll expresses the forecasting volatility also 
as a linear function of the past volatility and past returns. Eq.J22|) gives, however, the best linear forecast with the 
proved smallest forecasting error (|25p. The weight coefficients S T (0)' S ' allow to evaluate the magnitude and number 
of terms needed for the ARCH models to quantify future variance with sufficiently good precision. The ARCH models 
receive an additional support and more general framework through the BLF formula J22J). 

The accurate estimates of the future volatility are important for risk management and options pricing. The BLF 
formula (|23[) represents an interest as the proved most accurate estimate in the class of all linear functions of the past 
volatility and past returns. 
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